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ABSTRACT 

We present detailed model fits to observations of the disc around the Herbig Ae star HD 163296. This well-studied object has an age of ~4Myr, 
with evidence of a circumstellar disc extending out to ~ 540 AU. We use the radiation thermo-chemical disc code ProDiMo to model the gas and 
dust in the circumstellar disc of HD 163296, and attempt to determine the disc properties by fitting to observational line and continuum data. These 
include new Herschel/PACS observations obtained as part of the open-time key program GASPS (Gas in Protoplanetary Systems), consisting of a 
detection of the [Oi] 63 pm line and upper limits for several other far infrared lines. We complement this with continuum data and ground-based 
observations of the 12 CO 3-2, 2-1 and l3 CO J=l-0 line transitions, as well as the H 2 S( 1) transition. We explore the effects of stellar ultraviolet 
variability and dust settling on the line emission, and on the derived disc properties. Our fitting efforts lead to derived gas/dust ratios in the range 
9-100, depending on the assumptions made. We note that the line fluxes are sensitive in general to the degree of dust settling in the disc, with an 
increase in line flux for settled models. This is most pronounced in lines which are formed in the warm gas in the inner disc, but the low excitation 
molecular lines are also affected. This has serious implications for attempts to derive the disc gas mass from line observations. We derive fractional 
PAH abundances between 0.007 and 0.04 relative to ISM levels. Using a stellar and UV excess input spectrum based on a detailed analysis of 
observations, we find that the all observations are consistent with the previously assumed disc geometry. 

Key words. Stars: pre-main sequence; Circumstellar matter; Protoplanetary discs; Astrochemistry; Radiative transfer; Individual (HD 163296); 
Line: formation 


1. Introduction 

Protoplanetary discs set the initial conditions for planet forma- 
tion, and the processing of chemical species and dust grains as 
discs evolve may ultimately provide the seed for life. Indeed, 
conditions within circumstellar discs are conducive to a rich 
chemistry, with a large variety of conditions ranging from icy 
mantles on grain surfaces in the cold disc midplane to ionised gas 
bathed in stellar ultraviolet radiation at the disc surface. In order 
to attempt to probe the gas in discs it is necessary to observe gas 
emission lines from a variety of chemical species (molecules, 
atomic species, ions), and couple these observations with de- 
tailed thermochemical disc modelling. 

One important disc property for planet formation models is 
the amount of gas present as it evolves. The dust characteris- 
tics are generally better-understood than the gas, and the ISM 
gas/dust mass ratio of 100 is often adopted to give an estimate 
of the gas mass. However, gas masses derived from millimetre 
and sub-mm CO emission observations are typically lower than 
from dust observations with this assumption, sometimes by a 
factor ~ 100 dThi et al.[l2001i) . This is variously attributed to gas 
depletion due to pho toevaporation or planet formation, and CO 
freeze -out on grains (IZuckerman et al.Ul995l:lKaim) & Ber toldi. 
2000). These low level rotational CO transitions are also strongly 
dependent on the disc size, and often optically thick under cer- 
tain disc c onditions, limiting their ability to probe the dee per 
disc layers dPanic & Hogerheiidell2009llHughes et al.ll2008l) . In 
order to better constrain the disc gas mass it is necessary to ob- 
serve additional gas species, such as those observable in the far- 
infrared. The FIR range is a crucial observing window which 
can help to resolve some of these ambiguities. The fine struc- 


ture lines of atomic oxygen and ionised carbon probe the gas in 
the warm disc surface layers, and as products of CO photodis- 
sociation provide a good test for this proposed explanation of 
CO under-abundance. In addition, the high level rotational tran- 
sitions of CO arise from warm gas in the inner disc, probing 
different disc regions to the low level mm CO lines. 

In addition to the discrepancy between gas mass estimates 
from millimetre continuum and CO observations, there is uncer- 
tainty regarding the derived disc outer radii from CO and con- 
tinuum observations. The outer radius suggested by model fits 
to dust emission is often found to be smaller than that of the gas 
discyjn the case of the Herbig Ae star AB Aurigae. lPietu et aTI 
( 2005 ) suggested that a change in dust grain properties, leading 
to reduced opacity in the outer disc, could explain this discrep- 
ancy. For HD 1 63296, [Isella e t ail (12007 1) attributed the appar- 
ent difference in disc radii to a sharp drop in surface density, 
opacity or dust-to-gas ratio beyond 200 AU, while [Hughes et ali 
( 2008 ) proposed that the apparently conflicting dust and gas ob- 
servations could be reconciled using a different treatment of the 
disc surface density at the outer edge, motivated by similarity 
solutions for the time evolution of accreting discs. This study 
attributed smaller derived disc radii from dust observations to 
detector sensitivity thresholds, arguing in favour of larger dust 
discs, albeit with an exponentially-tapered density at the outer 
edge. A steepening of the density profile in the outer disc has 
also been observed in the case of DM Tau, LkCa 15 and MWC 
480 by Pietu et afl (120071) . 

The Herschel open-time key program GASPS (Gas in 
Protoplanetary Systems) aims to characterise the gas in discs at 
each stage of their evolution, by observing the far-IR lines of a 
wide range of discs, from young gas-rich accretion discs to gas- 


1 


I. Tilling et al.: Gas modelling in the disc of HD 163296 


Table 1 . Observed line data with Herschel/PACS. Detections are listed as Fi+cr) whereas non-detections quote < 3cr. There is an 
additional absolute flux calibration error of 30%. In the case of the [Cn] 1 58 gm line, strong contamination from the surrounding 
background prevented us from estimating a meaningful upper limit. 


Species 

A [p/m] 

v [GHz] 

Line Flux [lO-^Wm^] 

Continuum (RMS) [Jy] 

OI 

63.18 

4745.05 

193.1 ±5.8 

16.46 ± 0.07 

OI 

145.52 

2060.15 

< 8.5 

22.12 ± 0.03 

CII 

157.74 

1900.55 

- 

23.90 ± 0.05 

p-H 2 o 

89.99 

3331.40 

<9.4 

19.14 ± 0.04 

o-H 2 0 

179.52 

1669.97 

< 14.5 

22.28 ± 0.08 

0 -H 2 O 

180.42 

1661.64 

< 16.2 

21.58 ± 0.09 

o-H 2 0 

78.74 

3810.01 

< 15.0 

21.00 ± 0.10 

OH 

79.11 

3792.19 

< 17.0 

20.70 ± 0.10 

OH 

79.18 

3788.84 

< 17.0 

20.60 ± 0.10 

CO J=36-35 

72.85 

4115.20 

< 11.6 

17.30 ± 0.03 

CO J=33-32 

79.36 

3777.63 

< 22.8 

18.12 ± 0.07 

CO J=29-28 

90.16 

3325.12 

< 11.1 

19.10 ± 0.05 

CO J= 1 8- 1 7 

144.78 

2070.68 

< 13.1 

22.05 ± 0.04 


poor debris discs dMathews et all 1 2010 ). This paper represents 
the second in-depth study of a Herbig disc as part of this pro- 
gram, following after a stu dy focusing on the young Herbig Ae 
star HD 169142 dMeeus et al.Ll2010 ). 

HD 163296 is an isolated Herbig Ae star with spectral type 
Al Ve, mass ~ 2.5 M 0 and stellar luminosit y ~ 3 8 L 0 (this work), 
situated at a distance of 1 1 8.6^J^ ^pc (Ivan Leeuweni . |2007|) . 
Scattered light and millimetre continuum observations indicate 
the presence of an in clined circumstellar d i sc extending out to 
a rad i us of 540AU ( Mannings & Sargentt 119971: iGradv et akl 
2000 : Wisniewski et al. . l2008l) . In addition, there is evidence 
of an asymmetric outflow perpendicular to the disc, with a 
chain of six Herbig-Haro knots tracing its mass-loss history 
dDevine et al., 2000 ). The derived disc dust mass is in the range 


(5 - 17) X 1 (r 4 M,, dNatta etldlf2004tlTa nnirl<ula m et all [2008 a: 


iMannings & Sargent . 1997^ Isella et all 20071) . 

Recent near-infrared studies of the inner disc of HD 163296 
indicate an inner dust rim feature enclosing a bright emission re- 
gion extending inwards towards the star. This emission inwards 
of the dust rim has been attributed to an optically thin inner 
disc, although there is uncertainty regarding the size and com- 
position of such a disc. Derived radii for the dust rim lie in the 
range (0.2 - 0.55) AU ( Renard et all 1201 Ot lEisner et all 2009; 
Tanni rkulam et all l2008bl iBenistv et all 20101: iMonnieret al 


2006). An increase in opacity at the inner rim due to subli- 


mation of grains could cause the rim to puff up, and it has 
been suggested tha t this leads to time-variable self-shadowin g 
of the outer disc (ISitko et all 120081: IWisniewski et all 120081) , 
Mid-infrared imaging of warm dust in the disc surface layers 
seems to indicate little flaring dDoucet et all ’2006 ), and scat- 
tered light brightness profiles are consistent with a non-flared 
outer disc (IWisniewski et alil2008 ). Fits to the ISO-SWS spectra 
of HAeBe discs led lMeeus et al.1 ( 20011) to classify HD 163296 
as one of their group II objects, in which an optically thick inner 
disc shadows the outer disc from stellar radiation, resulting in a 
non-flared geometry. 

In this study we aim to fit a wide variety of observed emis- 
sion lines and the dust spectral energy distribution (SED) simul- 
taneously with a disc with this observed geometry. 


2. Herschel/PACS Observations 

We have obtained observations of HD 163296 in the far- 
infrared using Herschel/PACS, consisting of spectroscopic line 
observations and photometry. These include a detection of the 



Fig. 1 . Observed [Oi]63pun emission line with Herschel/PACS. 


[Oi] 63 pun fine structure line, and upper limits for a number 
of other atomic, ionic and molecular far-IR lines. We observed 
HD 163296 with the PACS spectrometer in line spectroscopy 
mode (obsid 1342192161, 1252 sec, 3 repetitions) and range 
scan mode (obsid 1342192160, 5141 sec, 3 repetitions), while 
chopping and nodding. 

The spectra were reduced using HIPE version 7.0, with 
the standard pipeline scripts for a chopped line scan of a 
point source. The reduction steps include division by the 
spectral response function and MedianOffsetCorrection, as 
the new flat fielding task is not yet robust for short range 
scans such as the present observations (PACS instrument 
team, private communication). For this paper, we only make 
use of the data contained in the central spaxel, to which 
we apply a wavelength-dependent correction factor for the 
loss of flux due to diffraction losses. With these additional 
steps, the absolute flux calibration accuracy is 30% (see 
PACS Spectroscopy Performance and Calibration Document, 
http://herschel.esac.esa.int/AOTsReleaseStatus.shtml). 

The line fluxes and upper limits are given with continuum 
data in Table Q] The PACS photometry data are given in Table [2] 
Comparison of the azimuthally averaged radial profiles with the 
observed point spread function (PSF), using Vesta, show that the 
disc is unresolved at both 70p/m and 160pm. The spatial extent of 
the observed [Oi] 63 pun emission is illustrated in Fig-El where it 
can be seen that the emission is dominated by the central 9. 4x9. 4 
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Table 2. Photometric fluxes observed with Herschel/PACS. 
These observations (AORS 1342228401 and 1342228401) are 
done with the standard miniscan. 


,1 [pm] 

70 

160 

Flux [Jy] 

18.5 + 3.3 

20.2+4.0 


arcsecond spaxel. This corresponds to a maximum radial extent 
of 560 AU for the emission. 

The pipeline reduction of the [Cn] 158 pm range shows this 
line in absorption, which is unexpected. Therefore, we also 
looked at the unchopped spectra in the different chop positions. It 
is clear that the region around HD 163296 is full of [Cn] 158 pm 
emission, and that this causes the emission feature at the position 
of the star to cancel out (and even become an absorption feature). 
This is likely due to cloud material along the line-of-sight, since 
HD 163296 is located close to the Galactic plane, and is consis- 
tent with the presence of CO emission at radial velocities offset 
from this object, but unresolvable with PACS dThi et al.Li200 1). 
We also inspected all the other lines for their unchopped appear- 
ance, but did not find contamination through the offset positions 
in those cases. 


3. Properties of HD 163296 

The study of the disc of HD 163296 is directly linked to the 
knowledge of the absolute parameters of the star itself and the 
behaviour of the spectral energy distribution (SED) in the ultra- 
violet (UV), optical and near infrared (near-IR), which is mostly 
dominated by stellar radiation; in particular, the energy emitted 
at wavelengths shorter than 2500 A is a fundamental piece of in- 
formation used to model the disc, due to its impact on the photo- 
chemistry and gas heating dWoitke et all 1 2010 :), This section de- 
scribes the stellar properties and studies the spectral energy dis- 
tribution at the range mentioned above. 

3.1. Stellar parameters 

The determination of the stellar parameters of HD 163296 was 
done in an extensive work devoted to the stu dy of a sample 
of Herbig Ae/Be stars dMontesinos et all l2009h . Details of the 
methodology followed and the observations used can be found 
in that paper. We outline in this section the relevant steps of the 
process. 

An iterative distance-independent algorithm based on the 
analysis of the SED and optical high-resolution spectra (for the 
estimation of the effective temperature and metallicity), and mid- 
resolution spectra of the Balmer lines H/j, Hy and Ho (for the 
estimation of the gravity), was applied. 

The knowledge of the effective temperature, gravity and 
metal abundance allows us to place the star in a logg, - log T e g 
diagram and superimpose the appropriate set of evolutionary 
tracks and isochrones for that specific metallicity ( lYi et al.L 
1200 ll) . This gives directly -or after a simple interpolation be- 
tween the tracks and isochrones enclosing the position of the 
star- the stellar mass and the age. Since there is a one-to-one 
correspondence between a pair (7' e fr, g*) on a given track and a 
pair (Teff, L*/L Q ), the stellar luminosity can also be estimated. 

Table [2 gives the stellar parameters of HD 163296 after a 
slight refinement of the results presented in iMontesinos et alJ 
(2003). 


Table 3. Parameters for HD 163296. 


Temperature 

9250+ 150 K 

Gravity (logg.) 

4.07 ± 0.09 

Metallicity ([Fe/H]) 

+0.20 + 0.10 

Mass 

2.47 ±0.10 M q 

Luminosity 

37.7 + 7.0 L 0 

Age 

4.2 ± 0.4 Myr 


Table 4. Photometry for HD 163296. 


Band 

A (pm) 

Magnitude and error 

Flux and error (Jy) 

U 

0.360 

6.95 ± 0.03 

3.00 ± 0.08 

B 

0.440 

6.92 ± 0.04 

7.27 ± 0.27 

V 

0.550 

6.86 ± 0.03 

6.56 + 0.18 

R 

0.640 

6.73 ± 0.05 

6.26 ± 0.29 

I 

0.790 

6.67 + 0.10 

5.48 + 0.51 

J 

1.215 

6.17 + 0.05 

5.55 ± 0.26 

H 

1.654 

5.49 ± 0.05 

6.69 + 0.31 

K 

2.179 

4.71+0.05 

8.56 + 0.39 


3.2. The SED: UV, optical and near-IR observations 

Simultaneous optical UBVRI and near-infrared JHK photome- 
try obtained as part of the E XPORT collaboration (lEiroa et all 
2001; Oudmaiier et al., 2001), were used to build the optical and 
near-IR photospheric spectral energy distribution of HD 163296. 
Table [4] shows the wavelengths, magnitudes and corresponding 
fluxes with their uncertainties. The calibration from magnitudes 
to flux es has been done us ing the zero points given by iBesselll 
d!979t) for the optical and ICoxl d 1999b for the near-IR magni- 
tudes. 

Ultraviolet IUE spectra of HD 163296 were extracted from 
the INES archivtJJ From the collection of 68 SW and LW (for 
“short” and “long wavelength”) spectra in low resolution, ob- 
tained through the large aperture of the spectrograph, all those 
that looked clean (few or none bad pixels) and with the cor- 
rect exposure classification codes were selected. The spectra 
(SW+LW) cover the range 1250-3000 A. 

A high-resolution, far-UV spectrum of HD 163296 was con- 
structed by combining data taken with HST STIS and FUSE. 
The STIS spectra, covering 1150 A to 3000 A, were obtained 
from the HST STIS Echelle Spectral Catalog of Stars (StarCafl). 
The FUSE spectrum of HD 163296 is extremely noisy below 
970 A; therefore, we used only the portion between that wave- 
length and 1150 A in our high-resolution composite UV spec- 
trum for HD 163296. FUSE spectra contain strong terrestrial 
airglow emission lines, due to the large size of the observing 
aperture. We removed these lines from the HD 163296 spectrum 
before stitching it to the STIS spectrum. 

A comparison of the IUE and FUSE+STIS spectra shows 
that they differ in the level of their continuum. The IUE spec- 
tra with the largest flux levels are a factor ~ 1 .6 higher than the 
FUSE+STIS spectrum, therefore, regarding the modelling, the 
ultraviolet contribution of the star will be represented by two dif- 
ferent UV input spectra, namely a “low UV” state, correspond- 
ing to the FUSE+STIS spectrum, and a “high UV” state, corre- 
sponding to the average of the highest IUE spectrr0 (see section 

ED . 

1 http://sdc.cab.inta-csic.es/ines/ 

2 http://casa.colorado.edu/~ayres/StarCAT/ 

3 The average IUE spectrum was build with the following individual 
spectra: SWP 28805, SWP 28811, SWP 28813 and LWR 02065. 
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RA (J2000) 


Fig. 2. [Oi |63//m spectra observed by Herschel/PACS towards HD 163296. Each spectrum corresponds to a 9.4x9. 4 arcsec pixel 
centered at the given coordinates. Emission is dominated by the central pixel, corresponding to a maximum outer radius of ~ 560AU, 
assuming a distance of 1 18.6 parsecs to this object. 


In addition, a near-infrared spectrum of HD 163296 was 
obtained at the NASA Infrared Telescope Facility (IRTF) on 
Mauna Kea on 2010 July 8 with SpeX (iRavner et all 12003b . 
Ten individual exposures of the target, each lasting 32 s, were 
taken using the short-wavelength cross-dispersed (SXD) mode 
of SpeX. This mode yields spectra spanning the wavelength 
range 0.8-2. 4 /jm divided into six spectral orders. The slit width 
was set to O'.' 3, which yields a nominal resolving power of 2000 
for the SXD spectra. Observations of HD 156717, an A0 V star, 
used as a “telluric standard” to correct for absorption due to the 
Earth’s atmosphere and to flux calibrate the target spectra, were 
obtained immediately preceding the observations of HD 163296. 
The airmass difference between the observations of the object 
and the standard was 0.06. The seeing was estimated to be ~ 
073 - 075 at 2.2 gm during the observations and conditions were 
clear. The statistical S/N varies across the spectral range but is of 
the order of several hundred over the entire SXD spectrum. The 
SpeX data were reduced with Spextool dCushing et all 120041) 
and the telluric correction was p erformed using the procedures 
and software described by iVacca et aE 1 2003 1. 

Figure [3] shows the SED of HD 163296. The solid black line 
is a photospheric model computed for T e $ = 9250 K, logg, = 
4.07 and [Fe/H]=+0.20 f rom the GAIA grid of mo dels created 
with the PHOENIX code dBrott & Hauschildtll2005h . The model 


has been reddened with E(B-V)=+0.15 (Ry = 3.1); this will be 
discussed in section I3T41 

3.3. Spectral type from the ultraviolet and near-infrared 
spectra 

Despite the fact that the precise determination of the stellar pa- 
rameters has been done through a careful analysis of a variety 
of observations, as we have mentioned in section [3TTl an attempt 
to determine the spectral type of HD 163296 in a more quali- 
tative way has been done using the ultraviolet and near-infrared 
spectra. 

A comparison of the ultraviolet data with IUE spectra of stars 
classified as “spectral-type standards” was done . The spectra 
were taken from the “IUE Ultraviolet Spectral Atlas of Standard 
Stars’ FliiwjT et al.Lfl983l.ll992l) . It is interesting to note that be- 
low 2000 A, for a given spectral type, the ultraviolet spectra of 
the “spectral-type standards” themselves are very different, in 
particular for A0 V and Al V. If we consider only the region 
above 2000 A, the comparison suggests that HD 163296 is closer 
to A3 V (and possibly slightly later), however, if a mild extinc- 
tion is present, the UV spectra would be also close to that of Al 


4 http://www-int.stsci.edu/~jinger/iweb/proj/project.html 
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Wavelength (ium) 

Fig. 3. Spectral energy distribution for HD 163296. The green 
dots represent the fluxes corresponding to the UBVRIJHK pho- 
tometry. The size of the dots is of the order or larger than the 
uncertainties. The IUE average spectrum and the FUSE+STIS 
spectrum are plotted in blue and red, respectively; the latter 
has been slightly smoothed to reduce the noise. The SpeX/IRTF 
spectrum is plotted in magenta. The black solid line is a photo- 
spheric model computed for the specific stellar parameters given 
in Table reddened with E(B-V)=+0.15 and normalized at the 
flux in V. See text for details. 


V stars. A larger extinction would push the spectral type deter- 
mination towards earlier spectral types. 

A similar exercise of spectral typing was done by a compari- 
son of the SpeX spectrum in t he YJ bands (0.85 - 1.35 pm) with 
Kurucz synthetic models (see lMontesinos et all (120091 ) for refer- 
ences regarding the synthesis codes and models used). The spec- 
trum is complex, the range of r e ff’s explored with the Kurucz 
models (8750-9250 K, corresponding to A3V and A1V) does 
not produce significant Paschen line variations. The Paschen line 
wings of HD 163296 are narrower than in any of the models, 
which would suggest temperatures of 9250 K or higher. The 
metallic lines are faint, their typical intensity being lower than 
5% that of the continuum. Some lines are deeper than in any 
photospheric spectrum, which suggests the presence of some 
circumstellar absorption. Some lines are much fainter, which 
suggests an spectral type A TV or earlier, or dust emission veil- 
ing. Good agreement is rarely found. Although we fail to find a 
unique combination of T e g and log g that provides an adequate 
fit to both the UV and NIR spectra of HD 163296, all data are 
consistent with an A1-A3 star, seen through mild foreground ex- 
tinction. 


Table 5. Fixed disc model parameters. 


Quantity 

Symbol 

Value 

stellar mass 

AU 

2.47 M e 

effective temperature 

F ft 

9250 K 

stellar luminosity 

U 

37.7 L e 

inner disc radius 

Kin 

0.45 AU 

dust material mass density 

Pgr 

3.36 gem - -* 

strength of incident ISM UV 

V ISM 

/V 

1 

cosmic ray Hi ionisation rate 

<fcR 

1.7 x 10- 17 s' 1 

turbulent Doppler width 

^turb 

0.15kms 1 

a viscosity parameter 

a 

0 

disc inclination 

i 

50° 

distance 

d 

118.6 pc 


R\r is not known. The observations available do not allow us to 
disentangle completely this problem. 

The strategy followed to assign a value of E(B-V) to this 
star has four steps: 1) introduce extinctions from E(B-V)=0.0 to 

0. 2 to the photospheric model, 2) normalize the reddened model 
photosphere to the flux at V, 3) deredden the normalized model 
and estimate the integrated stellar flux, Fp note that F* is the flux 
that would be observed in the complete absence of extinction. 
Obviously, different extinctions imply different values of F». A 
value of/? v = 3.1 has been used. 

The fourth step is based on the following argument: since 
the stellar luminosity L»/L 0 is known through the distance- 
independent method outlined in section 13.11 an estimation of 
the distance, d, can be obtained simply from the expression 
L» = 4nd 2 F », where it is assumed that the whole star is visible, 

1. e. and no part of the stellar flux is completely blocked by the 
disc of othe r circumstellar ma t erial ( see a complete discussion of 
this issue in lMontesinos et all (120091) '). The value E(B-V)=+0. 15 
allows us to recover a distance of 128 pc, which matches the 
Hipparcos distance to within the uncertainties, therefore it has 
been adopted as the optimum extinction for consistency with the 
whole set of parameters. We would like to point out that that 
value is a simple parameterization of a very complex problem 
and that a deeper study, outside the scope of this paper, would be 
needed to study all its peculiarities. 

In Fig. [3] the photospheric model reddened with E(B- 
V)=+0.15 has been plotted. It goes through the two ultraviolet 
spectra but falls below both of them at wavelengths shorter than 
~ 1260 A, and it also underestimates the flux at U. Note that 
the model is purely photospheric, therefore the excess fluxes in- 
ferred from the observed data when compared with the model 
might be due to emission from the accretion shock, i.e. with 
a non-photospheric origin. The observed SED seems to depart 
from the photospheric SED at around ~ 1 pm, where the contri- 
bution from the disc starts to become noticeable. 


3.4. The extinction 

The determination of the extinction is an intricate problem. The 
classical method of comparing the observed colours with those 
of standard stars of similar properties does not lead to any con- 
clusive result, given the slight variability of the star and the pho- 
tometric (and spectral type) uncertainties. For normal stars, the 
UV part of the SED can be used to estimate this parameter, but 
in the case of HD 163296 the ultraviolet spectrum is variable 
and has contributions both from the photosphere and the accre- 
tion shock that are difficult to quantify. The extinction seems to 
be mostly circumstellar, therefore the exact value of the constant 


3.5. Observational constraints 


In addition to the new Herschel data we have further con- 
strained our disc models using SMA interferometric observa- 
tions of the .1=3-2 transition in 12 CO and PBI interferometry of 
the I2 CO 2-1 and I3 CO 1-0 transitions ( Isella et al l l2007t) . an 
upper limit for the S(l) tra nsition in H 2 obtained by VLT/VISIR 
( Mart in-Zai'di et all 20101) . the ISO-SWS spectrum in the in- 
frared, and a wealth of photometric data, including SCUBA pho- 
tometry in the sub-mm dSandell et alll20l~fi) . See Fig.[4]for a full 
list of the photometry sources. 

There is a wealth of observational data available for HD 
163296 which constrains to varying extents our choice of model 
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parameters. As well as fitting to the data mentioned in the pre- 
vious paragraph, we have taken care not to contradict current 
understanding of the disc geometry. There is at present some un- 
certainty regarding the radial extent of the disc. Scattered light 
imaging indicates an outer edge for the dusty disc of ~ 540 AU, 
consistent with that derived for the gas disc from millimetre CO 
emission bv ilsella et alJ ( 2007 1. However, the resolved mm con- 
tinuum em ission seems to indicate a smaller oute r disc radius 
~ 200 AU (llsella et al.ll2007h . lHughes et alJ(l2008l) propose that 
the millimetre continuum and CO 3-2 emission can be fit simul- 
taneously by a disc with an exponentially-tapered outer edge. 
For the purposes of our modelling, we adopt the surface den- 
sity profile derived by [Hughes et al . (2008.) for this object, and 
attempt to fit simultaneously the remaining observational data 
(SED, line fluxes, CO line profiles) by varying the remaining 
disc parameters. As an additional test, we compute a model 
with a power-law density profile, varying in addition the column 
density power index to fit the same data (SED, line fluxes, CO 
line profiles). For the purposes of this power-law model the disc 
outer edge is fixed at the scattered light radius of 540 AU. These 
adopted disc outer radii are consistent with our far-infrared pho- 
tometry measurements, in which the disc is unresolved, corre- 
sponding to the emission at these wavelengths being dominated 
by the inner ~ 200-300 AU of the disc. 

The disc is also resolved at the inner edge, and while 
there is some uncerta i nty re g arding the precise st r ucture in this 


region dRenard et all 2010t Benistv et al.L 120101: lEisner et all 
2009|; iTannirkulam et all l20 08h). we adopt the median value 
of 0.45AU for the inner disc radius. The disc is assumed to 
have constant flaring with radius, i.e. the assumed inner rim 
structure does not cause shadowing. The reference scale heights 
referred to in our results (see Table 0 are defined at the in- 


constrained by imaging at various wavelengths (Benistv et al. 

2010; 

Wassell et al. 

2006; Isella et aljj2007; Tannirkulam et al.. 

2008a 

; Gradv et al. 

, 2000), and is fixed at 50° throughout the 


modelling in this paper. See Table [5]for a full list of the parame- 
ters which are fixed during the model-fitting process. 

While we have not used the spatial CO emission maps as a 
constraint for our modelling efforts, we do use the extracted line 
profiles for the three CO lines (llsella et all 120071) as a constraint 
for the radial origin of the emission. The interferometry allows 
us to isolate the emission from the disc, eliminating contami- 
nation from cloud material. The bulk of the modelling uses a 
surface density profile derived from fitting to the spatial CO 3-2 
emission (iHughes et all l2008h . 


4. Modelling 

The disc modelling was carrie d out using the radiation thermo- 
chemical disc code ProDiMo (IWoitke et all l2009t iKamn et all 
l2010l> . The code solves in turn the radiative transfer problem, 
chemical network, and gas heating and cooling balance. Finally, 
the level populations are calculated, followed by line transfer 
calculations to give the predicted line emission. 

The frequency-dependent 2D radiative transfer solver calcu- 
lates the dust temperature structure and internal continuous radi- 
ation field for a given disc and stellar spectrum. For the purposes 
of this paper it was used to determine an SED-fitting dust model 
for HD 163296, by varying the total dust mass, grain size pa- 
rameters and the spatial distribution of dust in the disc, including 
dust settling as a function of grain size (see Section l5~i~l >. 

Once the dust temperature structure 7jj ust (r, z) and internal 
radiation field J v {r,z ) have been computed, ProDiMo solves the 


chemical network assuming kinetic equilibrium (such that the 
net concentration of each species does not change over time), 
and the gas energy balance. This gives the chemical compo- 
sition of the disc, the gas temperature T gas (r,z) and the gas 
emission lines, and depends on the gas-to-dust ratio and PAH 
abundance. The chemical network includes 973 reactions be- 
tween 76 gas phase and solid ice species comp osed of 10 ele- 
ments (IWoodall et al.Ll2007tlSchoier et all [20051) There is a de- 
tailed treatment of UV-photoreactions (see Kamp et al., 2010), 
as well as Hj formation on grain surfaces, vibrationally excited 
H* chemistry, and ice formation (adsorption, thermal desorp- 
tion, photo-desorption, and cosmic -ray desorption) for a limited 
number of ice species (see IWoitke et all 12009L for details). 

The level populations of the various atoms, molecules and 
ions are calculated using an escape probability method, and 


case those of CO (Flower, 200 

\ Jankowski & Szalewicz 

2005 

Wernli et al., 2006; Yang et al. 

2010), HiO (Barber et al. 

2006 

Dubernet & Grosjean, 2002F 

Faure et all 2007S; Green et al. 


19931 
Hell et al„ 


12007 


1998; Chambaud et al., 1980; Jaauetetal 


1992 

19771 


Launav & Roue 
Launav & Roue 


1977) and Cn ( Flower & L aunav . 
19771 IWilson & Belli 120021) . See 


Woitke et al.l (201 1 ) for a description of the line transfer calcu- 
lations, as well as recent improvements to the chemistry and gas 
heating/cooling balance in ProDiMo. 


4.1. Input spectrum 

The stellar input spectrum used for the modelling is a com- 
posite of the available observed UV data with a model pho- 
tosphere spectrum in the wavelength range for which detailed 
spectral data are not available. We have accounted for the ob- 
served ultraviolet variability of this object (see Section. 0 by 
running two sets of disc models with different UV input spec- 
tra, namely, one representing the “low UV” state and the other 
a “high UV” state, as described in Section lT2l These were pro- 
vided by the FUSE+ST1S and the IUE average spectra, respec- 
tively. At the upper wavelength boundary we switch in each 
case to the GAIA PHOENIX model photosphere computed with 
T eS = 9250 K, log g* = 4.07 and [Fe/H]=+0.20. The observa- 
tions were de-reddened using the Fitzpatrick parameterization 
(Fitzpatrick, 1999) with E(B-V)=+0.15 and R v = 3.1. 

The UV spectrum is of central importance to the gas mod- 
elling, both in the chemical network and in the gas heat- 
ing/cooling balance. This requires UV input for wavelengths of 
912A and above (the wavelength at which atomic hydrogen is 
ionised). The FUSE spectrum in this case only extends down to 
970A, and the IUE spectrum down to 1150A. In order to give 
a full spectrum for modelling purposes, the “high UV” input is 
extended down to 970A by scaling the STIS+FUSE (low UV) 
spectrum accordingly. In both cases a power-law fit to the data is 
employed for the remaining shortfall. The total photon particle 
fluxes are 60% greater for the high UV input spectrum than for 
the low UV state. 


4.2. Fitting procedure 

An evolutionary ^-minimisation strategy was employed to si- 
multaneously fit the observed SED, ISO-SWS spectrum and gas 
line observations (see Section [3T5l >. In the case of the line obser- 
vations the models were fitted to the observed fluxes, and the line 
widths in cases where profiles were available. The upper limits 
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were not used as a fitting constraint, but were checked against 
the predicted best-fit model fluxes for any contradictions. 

We have utilised a Monte-Carlo evolutionary method, vary- 
ing the 1 1 parameters listed in Table [6] until the minimum x 2 
value was reached (see IWoitke et al.l (1201 ll) for more informa- 
tion regarding the evolutionary strategy and method of x 1 cal- 
culation). We note that the best-fit parameters depend in general 
on the weighting scheme employed for the three groups of ob- 
servations, and in practice it was found that weighting slightly 
in favour of the line emission and ISO-SWS observations at the 
expense of the photometry gave the best fit overall. We also note 
the possibility that the converged parameter values correspond to 
a local x 2 minimum, and not a single global minimum, and we 
do not claim a unique fit to the observations. Indeed, it is clear 
from our modelling efforts that there exist a lot of degeneracies 
between the various parameters. However, each of the models 
discussed in the following sections are the result of several hun- 
dred generations of ^-minimisation. 

4.3. Dust properties 

The rich solid-state spectru m of HD 163296 was observed by 
Ivan den Ancker et al.l d2000t) . using the ISO-SWS and ISO-LWS 
spectrometers. This study found the dust in this object to consist 
of amorphous silicates, iron oxide, water ice and a small frac- 
tion of crystalline silicates, with the presence of large millimetre- 
sized grains indicated by the continuum temperature. 

We assume homogeneous and spherical dust grains (Mie 
theory), with a power-law size-distribution defined by a mini- 
mum radius a mm , maximum radius a max , and power-law index 
p. The grain composition is assumed to be constant throughout 
the disc^ and we adopt the dust species mixture determined by 
iBouwman et alJ (i2000l l for this object, averaging their fractional 
species abundances over the disc (see Table [7]). All grains have 
the same species composition regardless of their size, and the 
grain composition is not treated consistently with the physics 
and chemistry of grains. For instance, the assumed grain water 
ice fraction is constant throughout the disc, independent of the 
water ice abundances calculated in the solution of the chemical 
network. For the dust material mass density, we take the aver- 
age value for this mix of 3.36 gcm^ 3 . Optical constants for the 
various dust species were taken from measurements made by the 
studies listed in Tabled 

[Bouwman et alJ (120001) carried out a detailed study of this 
object, fitting to spectral data over a large wavelength range. We 
note however that the dust opacity law represents a potentially 
large source of uncertainty when deriving the disc parameters. 
For instance, the two olivine species FeMgSKTt and Mg 2 Si 04 , 
whilst having broadly similar spectra in the wavelength range 
observable by ISO-SWS, have a factor > 30 difference in ab- 
sorption opacity at ~ 1 micron. This will inevitably affect the 
derived disc parameters, although the difference in opacities at 
millimetre wavelengths is negligible. The models referred to 
in this study have extinction opacities in the range (10.4-12.8) 
cm 2 g _1 (dust) at 1mm. 

4.4. Gas/dust ratio 

The models considered have constant dust properties with ra- 
dius. However, in general the models do allow for differential 
dust grain scale heights as a function of grain size, to represent 
the major effects of dust settling. This leads to a local gas/dust 
ratio which varies with height in the disc, yet at any given ra- 



X [|Xm] 


Fig. 4. Spectral energy distribution for the best-fit fully mixed 
disc model (“preferred” model), obtained from a simultaneous 
fit to the observed SED, ISO-SWS spectrum and line emis- 
sion from HD 163296 (solid black line). Black dotted line in- 
dicates the SED for the model with power-law density profile, 
which is unable to fit the mm continuum image for HD 163296. 
Blue circles indicate (with increasing wavelength) simultaneous 
UBVRIJHK photometr y dEiro a et all |2(Xnl lOudmaiier et all 
1200 ll) . LM photometry dde Winter et all 1 200 ll). IRAS pho tome- 
try (12-100 mic), sub-m m photome t ry dMa nnings. 19941) . and 
millimetre photometry (llsella et all 120071) . Also marked are 
PACS photometric observations (red circles), PACS continua 
derived from the spectroscopic observations (black circles), 
SCUBA photometry (green circles) (fSandell et alll20ll ). scaled 
VLT/UVES spectrum (blue line, Martin-Zaidi, in preparation) 
and the ISO-SWS spectrum (green line). The red line shows 
the stellar+UV input spectrum. Downwards arrow denotes up- 
per limit. All fluxes were corrected for interstellar reddening 
using the Fitzpatrick parameterisation (iFitznatrickl. 11999!) with 
/?y = 3.1 and E(B-V)=0.15. 


dius the ratio of the vertical gas and dust total column densities 
is constant, and equal to the global gas/dust ratio. The degree 
of dust settling is determined by a simple parameterisation (see 
Section [ATI) . and models in which no dust settling is present are 
referred to as “fully mixed”. 


5. Results 

Table [6] gives the derived parameter values for four runs of the 
evolutionary scheme outlined in the previous section. The first 
three columns refer to models with a fixed exponentially-tapered 
density profile. The fourth column refers to a model with power- 
law density profile, Xoc in which the outer radius is fixed at 
540 AU but the power index, e, is allowed to vary. 

The surface density profile in the models with an exponential 
outer edge (columns 1-3) takes the form 


I.(R)ccR 7 exp 



( 1 ) 
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Table 6. Model parameters for four well-fitting disc models. The first three columns refer to models with exponentially-tapered 
density profiles, while model 4 is a simple power law. Models 1, 3 and 4 are for a “low UV” input spectrum while model 2 is for a 
“high UV” flux. Finally, models 1, 2, and 4 allow dust grains to settle towards the midplane while model 3 is fully mixed. Model 3 
represents our overall “preferred” model, which is later discussed in detail, and is marked by an asterisk. 




CD 

(2) 

(3)* 

(4) 

Quantity 

Symbol 

“low UV” 

“high UV” 

Fully mixed 

Power-law density profile 

disc mass 

^clisc 

2.2 x 1(T 2 M 0 

2.0 x 10~ 2 M e 

7.1 x 10 - 2 M 0 

1.1 x 10 - 2 M e 

column density power index 

6 

n/a 

n/a 

n/a 

0.085 

reference scale height 

Ho 

0.019 AU 

0.019 AU 

0.019 AU 

0.027 AU 

flaring power 

P 

1.068 

1.066 

1.066 

1.019 

gas-to-dust mass ratio 

ptPd 

21.8 

20.5 

101.1 

9.1 

minimum dust particle radius 

^min 

5.95 x 10~ 3 pm 

1.25 x 10 -3 pm 

9.63 x 10~ 3 pm 

1.24 x 10~ 2 pm 

maximum dust particle radius 

^max 

11 34 pm 

1405 pm 

2041 pm 

1015 pm 

dust size dist. power index 

p 

3.61 

3.57 

3.68 

3.75 

dust settling parameter 

S 

0.57 

0.39 

0 (fixed) 

0.56 

minimum dust settling radius 

^set 

1.09 pm 

0.44 pm 

n/a 

0.81 

PAH abundance relative to ISM 

/pah 

4.3 x 10~ 2 

3.5 x 10~ 2 

6.8 x 10~ 3 

2.2 x 10“ 2 

fit to observed photometry 

Aphot 

1.723 

2.061 

1.323 

1.279 

fit to ISO-SWS spectrum 

Also 

1.326 

1.378 

1.780 

0.827 

fit to observed line data 

Aline 

0.661 

0.741 

0.742 

0.607 



Fig. 5. 1 .3mm continuum emission maps for HD 163296. From left: observed map; emission map for preferred model (column 
3 in Table [6}; residuals for preferred model; emission map for power-law model (column 4 in Table [6}; residuals for power-law 
model. Residuals computed by subtracting the model intensity from the observed intensity. Contours spaced at 12mJy intervals, 
corresponding to [3,6,9,12,15,18,21,24,27,30,33,36,39,42,45]x a. 


Table 7. Dust grain composition. 


Dust Species 

Mass Fraction Optical constant ref. 

Amorphous FeMgSi0 4 

0.745 

Dorschner et al. (1995) 

Amorphous Carbon 

0.15 

laser et al. ( 1998) 

Crystalline Mg 2 Si0 4 

0.035 

Servoin & Piriou (1973) 

Water Ice 

0.035 

Warren & Brandt (2008) 

Iron Oxide 

0.02 

Hennins et al. (1995) 

Iron 

0.015 

Posch et al. (2003) 


where y is analagous to the power law index, and Rq is the 
scale length over which the disc surface density tapers exponen- 
tially. We adopt the previously derived values of y = 0.9 and 
Ro — 125 AU (iHughes et al.L 120081) . For the purposes of our mod- 
elling, the disc chemistry etc. was computed out to 850 AU, at 
which point the column density is negligible. 

The first, “low UV” column uses the FUSE+STIS UV spec- 
trum as input for the modelling, and the second “high UV” col- 
umn uses the average IUE spectrum, as described in Section l4~Tl 
Both runs produce settled discs, with variable scale heights for 
dust grains of different sizes. The third column in Table[6]gives 
the results for a run with well-mixed dust grains, i.e. no settling 
is introduced. This run uses the FUSE+STIS (low UV) spec- 
trum, as does the run with power-law surface density referred to 
in column 4. 


None of the models described in Table [6] represent a perfect 
fit to all the available data for HD 163296, but in all cases the 
models are able to fit almost all the observations. Considering 
the large number of data used for the modelling, and the range 
of parameters covered by the well-fitting models, is is clear that 
there is some degree of parameter degeneracy present. 


5. 1 . Continuum emission and spectral energy distribution 

The SED for the best-fit fully mixed disc model (column 3 in 
Table[6j is shown in Fig. [4] With the constraints outlined in pre- 
vious sections, our aim was to fit the observed SED with a sim- 
ple dust model, i.e. a continuous disc with constant flaring, and 
dust species composition constant throughout the disc. All of the 
models described in Table [6] require the presence of mm-sized 
grains in order to fit the observed SED. It was not possible to fit 
both the 10 micron silicate feature and the millimetre tail with 
a well-mixed disc. This can be seen in Tabled where the fully 
mixed model gives a worse fit to the ISO-SWS spectrum than for 
the models in which dust settling is present. Fig.[4]shows that the 
best-fit fully mixed model gives a_smaller than observed silicate 
emission feature. ISitko et al] d2008l) observe a spectral variabil- 
ity of ~ 10% in the wavelength range covered by the ISO-SWS 
data, and while we are unable to fit the observations to within 
this range with even a settled model (simultaneously with the 
overall SED and line data), the fit is considerably better than 
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Fig. 6. Plotted in black are the line profiles observed by 
llsella et afl (i2007t) for the I2 CO J=3-2 (upper panel), 12 CO J=2- 
1 (middle panel) and 13 CO J=l-0 (lower panel) transitions, with 
the corresponding profiles from the preferred disc model in red. 
This refers to a single simultaneous fit to the observed contin- 
uum and line data (column 3 in Tabled. Observed profiles are 
obtained by integrating over the whole disc. 


is possible with a well-mixed disc. This result is similar to that 
found for IM Lupi bv lPinte et al.l (i2008h . in which a settled disc 
was needed to fit the 1 0 pm silicate feature and SED millimetre 
tail simultaneously. 


On the other hand, the constraint of fitting to the ISO-SWS 
spectrum leads to models which produce too much emission in 
the near-infrared (J,H & K bands), giving an overall worse fit 
to the observed photometry. This effect is less pronounced in 
the best-fit fully mixed model (column 3 c.f. columns 1&2 in 
Table [6|. The flux overprediction in the settled models is likely 
mainly due to the smaller average grain size in the strongly il- 
luminated disc surface layer. The condition of radiative equilib- 
rium between grains means that this gives a higher average grain 
temperature at the disc surface, and the re-emitted light peaks 
at shorter wavelength than for the fully mixed models (in the 
near-IR as opposed to the mid-IR). This flux overprediction is 
still present to a smaller extent in the fully mixed model (see 
Fig. 0). This could be a consequence of our fixing the disc in- 
ner radius for modelling purposes, and a better fit in the near-IR 
would be possible if we had allowed the inner radius to vary. It 
is also probable that the structure of the inner rim is more com- 
plex than allowed by our parameterised model, with evidence 
for puffing-up of material, which would certainly affect the emis- 
sion at these wavelengths. There is also significant variability ob- 
served for this object in this wavelength region (Ide Winter et all 
12001 IlSitko et~aTll2008h . 

Dust settling enhances the silicate emission feature for an 
optically thick disc since it removes the large grains from the 
surface layers, and places them at smaller heights in the disc 
where they cannot be observed at mid-infrared wavelengths. The 
flat blackbody opacity of these larger grains is overwhelmed by 
the characteristic spectrally-varying opacity of smaller micron- 
size grains which remain at lower optical depth, dominating the 
observed emission. 

In cases such as this where a simple parameterised disc struc- 
ture is assumed, ProDiMo implements a simple recipe to ac- 
count for the major effects of vertical dust settling. We assume 
that the dust grains are distributed vertically with a scale height 
which decreases for large dust particles: 

H'(a,r)-H(r) ■ maxjl, a/a s }~ s ^ 2 (2) 

where H(r) is the gas s cale he ight , and s and a s are two free 
parameters (see lWoitke et al.ld2010b for details). 

Introducing two additional parameters to the modelling in- 
evitably improves the fit to the observed data, and indeed the 
combined fit to the photometry, ISO-SWS spectrum and line 
emission does improve overall as dust settling is introduced. 
The fact remains, however, that it is possible to obtain a good 
fit to almost all the available data with a well-mixed disc, with 
a gas/dust ratio almost equal to the canonical value of 100, and 
so we adopt this model, described in column 3 of Table [6] as 
our “preferred” model. However, the inability of a well-mixed 
disc such as this to fit the observed emission in the mid-infrared 
means that there remains compelling evidence for a disc exhibit- 
ing dust-settling, able to fit the line emission with a depleted 
gas/dust ratio (see column 1 in Tabled. The effect of dust set- 
tling on the line emission is further explored in Section [5.2.2l 

The model with power-law density profile (column 4 in 
Table [6)1 gives the best overall fit to the SED, ISO-SWS spectrum 
and line emission (see Fig. [4}, with quite different disc parame- 
ters to the three models with exponential outer edges (columns 
1-3). The disc requires a very flat density profile, e ~ 0.085, 
in order to fit these observations (not including spatial data). 
The predicted millimetre continuum emission from this model 
is plotted with the observed maps in Fig. [5] It is clear that the flat 
density profile leads to an emission deficit in the inner disc in 
comparison to observations, and so this cannot be considered to 
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be a good model for the disc of HD 163296. A better power-law 
model could be found by using the maps as a further constraint 
for the modelling, but that is beyond the scope of this paper. 
The fully mixed “preferred” model gives a better match to the 
observed millimetre emission, as should be expected since its 
density profile is derived from fitting to this data dHughes et al.L 
2008 ). The factors driving the power-law model to such a flat 
density profile will be discussed in Section 15721 

The power-law density profile model is, however, able to 
fit the non-spatially-resolved data with a smaller flaring in- 
dex (~ 1.02) than the exponentially-tapered discs (~ 1.07), 
more in keeping with the observational evidence for a disc 


that is not strongly flaring (S 

eeus et alj. 2001: Doucet et all 

2006; Wisniewski et al.l 2008). 

Meeus et al. (2001) classify HD 


163296 as a group II object, in which the inner disc shadows 
the outer dis c, as opposed to group I objects which have a flared 
geometrv. lMeiier etail (12008 ) find HD 163296 to lie on the tran- 


sition between flared and non-flared geometry, i.e. a flaring index 
close to 1 . It is clear that our “preferred” model is unable to ac- 
count for every observed property of the disc of HD 163296, and 
it is probable that a more complex model - with non-constant 
flaring and variable dust properties with radius - would be re- 
quired to simultaneously fit the full set of observational data 
available for this object. 

The models described in Table[6]have disc dust masses in the 
range (7 - 12) X 10~ 4 M o , which is within the range of masses 
(5 - 17) x 1 0 4 M 0 , found in the literature for HD 163296. This 
small spread in dust mass from our fitting efforts is unsurprising 
given the fixed disc size and constant grain composition. 


5.2. Gas properties 


The predicted line fluxes for the various models are given in 
Table 0 We have ignored the [Cn] 1 58 pm result during the 
model-fitting process, due to systematic uncertainties arising 
from strong emission at the offset positions. 

The observed line profiles for the 12 CO J=3-2, 12 CO J=2-l 
and 13 CO J=l-0 transitions (llsella et all 1 20071) are plotted with 
the preferred model profiles in Fig. [6] The observed profiles are 
double-peaked, consistent with a disc in Keplerian rotation. The 
CO 3-2 line is observed to be brighter than predicted by the 
model. This behaviour is repeated across the three exponentially- 
tapered models, but is less pronounced in the power-law model, 
whose flat surface density profile gives a larger CO 3-2 flux. 
Discs with flatter density profiles allow stellar radiation to pen- 
etrate deeper into the inner disc, resulting in a slight increase 
in temperature throughout the disc, including its outer regions. 
This leads to brighter CO 3-2 lines since this transition is opti- 
cally thick in all of our models, and is largely dependent on the 
disc outer radius (which is fixed) and the outer disc temperature 
at intermediate height (see Fig. [8]). The CO 2-1 line also fol- 
lows this behaviour, with the line flux ratio in our models staying 
roughly constant at a value of CO 3-2/1 -0 ~ 3, as expected for 
lines formed under optically thick LTE conditions dKamn et al.L 
120101) . This is slightly lower than the observed ratio of CO 3- 
2/2-1 = 4.3. We note that recent observations of this object give 
a peak intensity of 25 Jy f or the 3-2 lin e ( Hughes et a l.. 201 10, 
compared with 43 Jy in the llsella et al.l d2007h data and the 30 Jy 
predicted by our preferred model. Our predicted flux is within 


5 This CO 3-2 data also seems to indicate t he presence of turbu lence 

of order 0.3km/s in the disc of HD 163296 dHughes et all [201 1). We 

have not yet explored the effect of this possible stronger turbulence on 

our results. 


the calibration error margin for both observations, but fitting to 
the lower value would in general tend to lead to power-law mod- 
els with steeper density profiles, in contrast to the flat power- 
law profile we obtain from fitting to the brighter observation. 
Another factor driving the power-law models to flat density pro- 
files is the overprediction in near-IR continuum emission by our 
models. The inner radius is fixed in accordance with high resolu- 
tion imaging, and so the power-law models tend to reduce their 
near-IR emission by removing material from the inner disc, giv- 
ing a better fit to the photometry. 

PAHs are one of the main sources of gas heating in the disc, 
as UV photons cause them to emit excited electrons via the pho- 
toelectric effect, which thermalise in the gas. For the purposes 
of our modelling we consider a typical size of PAH molecule 
(circumcoronene, Nq = 54 carbon atoms and Ah = 18 hydro- 
gen atoms) and include PAH , PAH, PAH + , PAH 2 _and PAH 3+ 
as additional species in the chemical network (see lWoitke et al.l 
(1201 Oh for further details). The fractional PAH abundance, /pah, 
is defined relative to the standard ISM particle jibundance with 
respect to hydrogen nuclei, Xp™ =3 x 10 5 * 7 * (iTielensl .2008), so 
that the PAH abundance, e(PAH) =/pah -^pah 7 ^- S ma H /pah val- 
ues (0.007-0.04 relative to ISM abundances) are required in our 
models in order to fit the observed line data. This is consistent 
with the findings of Geers et a l.l ( 2006 ). who require PAH abun- 
dances /pah <0.1 in order to fit the observed PAH emission from 
Herbig discs. PAH emission was not observed in HD 163296 by 
lAcke et akl d2010h . an d we have used th e Monte Carlo radiative 
transfer code Mcfost dPinte et all [2006 ) to check that our results 
are consistent with this non-detection. A disc model with param- 
eters and PAH abundance identical to our preferred model gives 
an infrared spectrum in which no PAH emission is seen. This is 
due in part to the low flaring in the disc, / = 1 .066, and the large 
fraction of small dust grains in the disc, with significant opacity 
in the optical-UV range over which PAH molecules absorb. This 
results in the PAHs being “hidden” from direct stellar illumina- 
tion, absorbing a fraction ~ 1 0 4 of the energy absorbed by the 
disc. In our preferred model 0.3% of the total carbon mass is tied 
up in PAH molecules. 

The total hydrogen number density and gas temperature 
within the disc are plotted in Fig. |Tj and the spatial origin of 
the various emission lines is visualised in Fig. [8] These plots all 
refer to the “preferred” disc model, i.e. one in which no dust 
settling is present. The CO .1=3-2 and [Oi] 63 pm lines are op- 
tically thick throughout the disc, and cannot be used alone to 
trace the gas mass. Even the l3 CO line is optically thick through- 
out much of the disc, only becoming optically thin outside of 
~ 400 AU. Also, there is evidence that this line (and most oth- 
ers) can be affected by the degree of dust settling in the disc 
(see Section 15.2.21) . The [CII] 1 57.74/rm line is optically thin, 
but traces only the warm ionised gas in the disc surface. Clearly 
care must be taken when attempting to use individual emission 
lines as a tracer of gas mass. Note the contrast in disc radii at 
which line and continuum become optically thin in the case of 
CO 3-2, reflecting the conflicting derived radii for the gas and 
dust discs in HD 163296. The line optical depths in Fig. [0 were 
computed using an escape probability formalism, and do not take 
into account the effects of Keplerian shear in the disc. 

As well as optical depth effects, disc gas mass estimates 
from CO emission alone can be affected by the extent of CO 
freeze-out on grain surfaces. ProDiMo calculates the grain ad- 
sorption and desorption as part of the solution of the chemical 
network, and the spatial abundances of gas-phase and ice CO in 
our best-fit model are plotted in Fig. [9] The total mass of gas- 
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Table 8. Integrated line fluxes for the transitions observed by PACS, and additional CO and H 2 lines, as predicted by the models 
listed in Table[6] Observed fluxes for comparison. Fluxes in [10 13 W/m 2 ]. Errors are a quadratic sum of the calibration error (20% 
for the interferometric CO observations and 30% for the PACS observations) and the RMS continuum noise. 


Species 

A Lum] 

v [GHz] 

“low UV” model 

“high UV” model 

Fully mixed model 

Power-law 

Observed Flux 

OI 

63.18 

4745.05 

200.6 

222.0 

191.4 

170.3 

193.1(58.2) 

OI 

145.52 

2060.15 

6.29 

7.33 

5.39 

4.59 

<8.5 

CII 

157.74 

1900.55 

8.32 

8.40 

11.2 

9.98 

- 

P-H 2 O 322 — > 2ll 

89.99 

3331.40 

1.22 

1.68 

1.20 

0.353 

< 9.4 

o-H 2 0 2 12 — » loi 

179.52 

1669.97 

4.40 

4.19 

2.43 

3.94 

< 14.5 

0 -H 7 O 271 — > 2i7 

180.49 

1661.64 

0.847 

0.993 

0.618 

0.243 

< 16.2 

0 -H 7 O 4^? — > 3)7 

78.74 

3810.01 

2.33 

2.94 

2.20 

0.840 

< 15.0 

o-H 2 S(l) 

17.03 

17603.78 

1.15 

1.59 

1.15 

0.572 

<28 

OH 

79.11 

3792.19 

4.59 

6.31 

5.40 

2.71 

< 17.0 

OH 

79.18 

3788.84 

4.78 

6.54 

5.59 

2.82 

< 17.0 

CO J=36-35 

72.85 

4115.20 

0.335 

0.489 

0.423 

0.132 

< 11.6 

CO J=33-32 

79.36 

3777.63 

0.497 

0.704 

0.578 

0.197 

<22.8 

CO J=29-28 

90.16 

3325.12 

0.763 

1.03 

0.732 

0.311 

< 11.1 

CO J=18-17 

144.78 

2070.68 

1.91 

2.31 

1.49 

0.851 

< 13.1 

CO J=3-2 

866.96 

345.80 

1.26 

1.27 

1.22 

1.49 

1.65 (0.39) 

CO J=2-l 

1300.40 

230.54 

0.401 

0.401 

0.397 

0.464 

0.379 (0.118) 

13 COJ=1-0 

2720.41 

110.20 

0.0122 

0.0115 

0.0120 

0.0127 

0.0124(0.007) 



Fig. 7. The total hydrogen number density (left panel) and gas temperature (right panel) are plotted for a vertical cross-section 
through the preferred disc model (column 3 in Table[6]i. Dashed lines indicate contours of gas temperature and visual extinction Ay, 
as marked. 


phase CO in the disc is 1.51 x 1O~ 5 M 0 , and the CO ice mass is 
4.32 x 10~ 5 M g . We would expect a smaller fraction of CO to be 
frozen out on grain surfaces in discs with flatter density profiles, 
due to weaker dust shielding in the inner disc leading to warmer 
temperatures throughout the disc, and less CO freeze-out. 

The derived chemical abundances and gas properties in our 
preferred model have been checked by re-computing the chem- 
istry using a time-dependant solver, and comparing the results 
to those obtained with the assumption of kinetic equilibrium. 
We assume molecular cloud initial abundances, and after run- 
ning the solver for 4Myr there is no major departure from the 
equilibrium chemistry. The biggest change in line flux is a re- 
duction of 20% in the 180.4/rm water line, while the [Oi] 63 gm 
line decreases by 4% and the CO lines by < 1%. The assumption 
of a constant l3 CO/ 12 CO ratio is also valid since dust shielding 
dominates over CO self-shielding in the models, with negligible 
change to the results when CO self-shielding is swtiched off. We 


would therefore expect there to be no fractionation effect present 
in this disc. 

All of the disc models considered are passive discs, i.e. the 
viscous “a” heating parameter is set to zero. This is not strictly 
consistent with the presence of accreting gas in this object, since 
this implies some form of viscosity to transport angular momen- 
tum through the disc. However, viscous heating is likely to be 
unimportant in the case of Herbig Ae discs such as this, which 
are thought to be dominated by radiative heating by the central 
star fo’Alessio et all 119981 ). 

5.2.1 . Effect of UV variability 

We have investigated the effects of UV variability on the gas in 
the disc in two ways. Firstly, by computing a test model with 
parameters identical to our “preferred” model, but with the high 
UV state spectrum used as input. This allows us to isolate the 
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Fig. 8. Spatial origin of the various gas emission lines in the preferred model (column 3 in Table |6j. From top-left clockwise: CO 
J=3-2, CO J=18-17, [CII] 157.74/rm, 0 -H 2 S(l) (with red T g contours), [OI] 63.18^/m, 13 CO J=l-0. In each panel, the upper plot 
shows the line optical depth as a function of radius (blue line) and the continuum optical depth at the corresponding wavelength 
(black line). The middle plot shows the cumulative contribution to the total line flux with increasing radius. The lower plot shows 
the gas species number density, and the black line marks the cells that contribute the most to the line flux in their vertical column. 
The two vertical dashed lines indicate 15% and 85% of the radially cumulative face-on line flux respectively, i.e. 70% of the line 
fl2ix originates from within the two dashed lines. 
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Fig. 9. CO abundances in the preferred model (column 3 in Table [6]). Left panel: gas-phase CO abundances. White and blue dashed 
lines plot the gas and dust temperature contours respectively. Right panel: CO ice abundance. 


Table 9. Effects of UV variability. Comparison between the preferred model, and an identical model illuminated by an input 
spectrum typical of a “high UV” state for this object. Shown are total species masses, average overall CO temperature, average CO 
temperature for the warm gas regions in the inner disc, and line fluxes. Luv refers to the energy emitted in the wavelength range 
912-2500A. 



Best-fit model 

Test model (high UV) 

Tuv/L* 

0.097 

0.155 

M h 

2.67 x 10- 5 M o 

2.82 x 10~ 5 M o 

M c * 

8.67 x 10~ 8 M o 

9.42 x 10~ 8 M o 

Mqo 

1.54 x 10~ 5 M o 

1.66 x 10~ 5 M o 

(Tco) 

55.3 K 

55.0 K 

(T co )(z/R> 0.2, R< 20 AU) 

262.8 K 

332.6 K 

F [OIJ 63 fim 

1.91 x 10~ lfl W m~ 2 

2.49 x 10~ lb W nT 2 

Ff CII]58/mi 

1.12 x 10~ 17 W nr 2 

1.41 x 10~ 17 W nr 2 

Fco 3-2 

1.22 x 10~ I8 W m -2 

1.31 x 10~ 1S W m -2 

FcO 36-35 

4.23 x 10~ I9 W m~ 2 

6.22 x 10~ 19 W nT 2 


effect of increasing the UV intensity on the line emission, since 
all other model parameters remain the same. The results of this 
test are summarised in Table [9] The second approach consists of 
a separate evolutionary run of models using the high UV state as 
input, and has been covered in Section 0 This method allows us 
to assess how the disc model parameters change in order to fit the 
same line data with an increased UV intensity (see Table |6}. By 
examining firstly the physical effect of the UV on the gas in the 
disc, and secondly the effect of the UV on the model parameter 
fit, we can estimate the degree of uncertainty introduced by UV 
variability in this object. 

An increased UV strength affects the gas chemistry by pro- 
moting photochemical reactions, and the photodissociation of 
H 2 and CO. It also leads to increased desorption of ice species 
from grain surfaces in regions where the UV is able to penetrate. 
It also strongly affects the gas heating, via the photoelectric ef- 
fect in dust grains and PAHs. 

The “high UV” input spectrum represents a 60% increase in 
UV intensity over the low UV input. Table |9]summarises the ef- 
fects of this increase on the gas in the preferred model. There is 
a slight overall increase in the disc gas temperature, and an in- 
crease in the atomic hydrogen and C + masses, due to photoelec- 


tric effects and PAH heating, and increased photodissociation of 
gas molecules. The total gas-phase CO mass actually increases 
slightly, since the increased rate of photodissociation is balanced 
by an increase in the rate of ice desorption. The J=36-35 line in 
particular sees the biggest increase in flux, ~ 50%. By consider- 
ing only the gas in the warm layers of the inner disc, where the 
high J CO lines form, it can be seen that this increase in line flux 
is caused by a substantial increase in the gas temperature in this 
region. The [Oi] 63 /.im flux increases by ~ 30%, roughly equal 
to the calibration uncertainty margin, and the low J CO lines in- 
crease by ~ 5%. All non-detected lines stay within the observed 
upper limits. In summary, it would not be possible to distinguish 
between the high and low UV states from these disc observa- 
tions, since all fluxes are within the observational error margins. 
The general influence of the UV intensity on the disc gas is dis- 
cussed by IWoitke eFal . ( 1201 0). in the context of a large grid of 
models. 

The derived parameters for the evolutionary run using the 
“high UV” spectrum as input are listed in Table [6] They are 
largely similar to those derived from the low UV run (columns 
2 c.f. 1), as would be expected from the small fractional change 
in UV intensity. The slight decrease in gas mass and flaring will 
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tend to counteract any increase in line emission. The PAH frac- 
tional abundance decreases slightly as might be expected. The 
main change lies in the grain size distribution, with a reduction in 
the minimum grain size. This makes it harder for the UV to pen- 
etrate the disc, balancing the increase in UV intensity. This also 
produces even more continuum emission in the near-IR, wors- 
ening the SED fit still further in comparison to the “low UV” 
model. 

5.2.2. Effects of dust settling 

The effect of dust settling on various line fluxes can be seen 
in Fig. [TO] This shows the increase in line flux across the var- 
ious transitions as smaller and smaller grains are allowed to set- 
tle towards the midplane. The models are otherwise identical to 
the preferred (fully mixed) model, but with a settling parame- 
ter of 0.5. The models run from strongly settled discs through 
to entirely well-mixed, i.e. identical to the preferred model. For 
the purposes of this exercise we have computed some extra line 
fluxes in addition to those observed in HD 163296, namely the 
12 CO 6-5 and 1-0 transitions, the 13 CO 3-2 transition, and the 
0 -H 2 O 538.3//m transition. 

iKamp et ail (1201 ll) computed a larg e grid of disc mod els us- 
ing ProDiMo, in tandem with Mcfost (Pinte et al., 2006). This 
study noted a trend of increasing disc dust temperature with set- 
tling, due to a combination of reduced emissivity at long wave- 
lengths in the disc surface, and increased illumination of the dust 
as the stellar radiation is able to penetrate further into the disc. 
Our models of HD 163296 follow this same behaviour, albeit 
with an accompanying increase in gas temperature. The effect 
of dust settling on the gas and dust temperature contours in our 
models is demonstrated in Fig. |TT] where the chemical abun- 
dances of various chemical species are plotted. The chemical 
structure is seen to “follow” the settled dust grains towards the 
midplane, since the UV penetrates further into the disc, caus- 
ing the characteristic layering of the various gas species to move 
closer to the midplane. This is despite the gas scale heights be- 
ing identical in both models. In the settled models the warm 
gas in the disc resides at lower heights where the gas density 
is higher, leading to a general brightening of the emission lines. 
This seems to be a firm result for this object, with serious im- 
plications for the derivation of a precise disc gas mass from the 
emission lines. 

The brightening of the line fluxes seems to be a general re- 
sult for all the transitions considered. The exact behaviour of the 
various line fluxes with variable dust settling, and the extent to 
which they are affected, depends on the height and radial posi- 
tion in the disc from which they originate. This is well-illustrated 
by the CO lines in Fig. [TO] As the largest grains begin to settle, 
the 13 CO 1-0 line is the first to show an increase in flux. This is 
due to an increase in thermal desorption of CO ice from grain 
surfaces in the outer disc midplane, which this line traces (see 
Fig. [8J. This behaviour is echoed in the other 13 CO line, the .1=3- 
2 transition. In contrast, the 12 CO 3-2 line, which is also formed 
in the outer disc, shows a smaller increase in flux. This is be- 
cause this line is optically thick, and formed higher in the disc 
where it remains unaffected by the thermal desorption of ice in 
the midplane. The 12 CO 2-1 and 1-0 lines follow this same be- 
haviour (not plotted). As smaller grains begin to settle, the 13 CO 
lines increase less rapidly than for the lines formed in the warm 
gas higher in the disc. This is most pronounced in the high J CO 
lines, which form in the inner disc where the temperature and 
density gradients are more extreme, and the lines more sensitive 
to the downwards shift in chemical structure. The 12 CO 6-5 line 



Fig. 10. The fractional enhancement in line flux is plotted against 
the minimum grain size affected by dust settling, relative to the 
preferred (fully mixed) model (a s = a max ). All models have a 
settling parameter of 0.5. Top panel: CO lines; Middle panel: 
water lines; Lower panel: [Oi]63/rm, [Cn]158/zm, OH 79.11/rm 
and Hi S(l) lines. 
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forms at intermediate distance in the disc, and this is reflected in 
the level of flux enhancement (Fig.flOl top panel). 

The behaviour of the water lines (middle panel in Fig. ITOt is 
less straightforward. The 0 -H 2 O 179.5/rm line shows a stronger 
increase in flux than both the 1 80.5/on line, which is at higher 
excitation, and the 538.3/rm line, which is at lower excitation. 
The P-H 2 O 89.99/tm and 0 -H 2 O 78.7/tm lines are at higher ex- 
citation, forming in the inner disc, and these also show a large 
increase in flux as the dust grains settle. It is likely that the rel- 
ative behaviour of the various water lines is a complex function 
of the changing temperature structure and H 2 O chemical abun- 
dance structure in the disc, affecting the excitation of the various 
levels. 

The largest increase in line flux occurs in the H 2 S(l) tran- 
sition (Fig. [TOj lower panel). The emission region for this line 
is seen to follow the 300K gas temperature contour (see Fig. ED, 
which intersects a greater molecular hydrogen mass in the set- 
tled models. This leads to a dramatic enhancement in flux in 
the strongly settled models, a factor of ~ 30 increase on the 
fully mixed model. The [Oi] 63 fim line also starts to increase as 
the smallest grains are allowed to settle, coinciding with an in- 
crease in atomic gas, but the enhancement is less extreme than 
since this line is formed further out in the disc (see Fig. [8]). The 
OH 79.11 gm line is formed at intermediate distance between H 2 
S(l) and [Oi] 63 gm, and its behaviour with increasing dust set- 
tling reflects this. The [Cn] 158//m line is less strongly affected 
by settling. This line is formed in a thin layer at the disc sur- 
face, with emission dominated by the outer disc, and so it is less 
sensitive to changes in the internal disc temperature structure. 

The settling of dust grains might be expected to drive line 
formation conditions closer to LTE, with the various species re- 
siding lower in the disc, in regions of higher density, with more 
frequent collisions between particles. However, there is no evi- 
dence for this in our models. The majority of lines show only mi- 
nor departures from LTE, in both settled and well-mixed models. 
The largest departure from the line fluxes assuming LTE level 
populations occurs in the water lines, but the LTE line fluxes in- 
crease at roughly the same rate as those computed using escape 
probability, so the ratio Fnlte/Flte remains roughly constant. 
The critical density for the water lines, n CI j t ~ 10 8 - 10 10 cirr 3 . 
These densities are reached only in the inner disc in our models, 
with each of the lines arising in regions with n < n cnl , even in 
the settled models. The water line which is closest to LTE in our 
models is the 538.3/rm line, which has the lowest level of ex- 
citation and the smallest critical density. In general, the critical 
densities decrease rapidly at large optical depths, driving most 
of the transitions towards LTE. 

These results are in contrast to the findings of tlonkheid et alJ 
(120071) . who find a trend of decreasing gas temperature and line 
emission as the dust grains are allowed to settle towards the mid- 
plane. However, a direct comparison is not appropriate, since 
their models assumed a simultaneous reduction in the dust/gas 
ratio and PAH abundances with dust settling. The reduction 
in gas temperature can largely be attributed then to the drop 
in photoelectric heating from dust grains and PAH molecules. 
Llonkheid et al.l ( 2007 ) also note that the dust temperatures in 
their models are high enough to prevent CO freeze-out, whereas 
considerable freeze-out is present here. iMeiierink et al.l ( j2009i) 
also predict an increase in line emission with dust settling, from 
models in which the effect of settling is simulated by simply in- 
creasing the global gas/dust ratio, as opposed to considering the 
vertical distribution of grains in the disc. 

Our models indicate that the 13 CO 1-0 line flux can change 
by a factor ~ 2 in strongly settled models. This is in a high mass 


disc for which the line is optically thick throughout most of the 
disc, and so this settling flux enhancement might be expected to 
be even stronger in a lower mass disc. In any case, variable grain 
freeze-out would seem to limit the ability of this line to trace the 
total disc gas mass. In general, the effect of dust settling on the 
vertical thermal structure of the gas in discs is seen to introduce 
new degeneracies when attempting to fit the disc parameters to 
the observed line emission. While our preferred model fits the 
observed data with a well-mixed disc with the canonical gas/dust 
ratio, we cannot rule out the possibility of a gas-depleted disc 
in which dust settling gives an enhancement in the various line 
fluxes. Indeed, such a disc with gas/dust ~ 20 allows a better 
fit to the line data and the observed 10 micron silicate emission 
(column 1 in Table [6]). 

5.3. Effect of X-rays 

It is currently unclear how important a role X-rays play in de- 
termining the gas chemistry and temperature structure in discs. 
lAresu et al. (1201 lh studied the effects of stellar X-ray emission 
on models of T Tauri discs computed with ProDiMo. This study 
found that while Coulomb heating by X-rays introduced an ex- 
tended hot gas surface layer to the discs, the [Oi] and [Cn] 
fine structure lines were only affected for X-ray luminosities 
Lx > 10 30 erg s 1 . One w ould expect this to hold f or the warmer 
gas in Herbig discs, and Gunther & Schmitt (2009) derive an X- 
ray luminosity of Lx = 10 29 6 erg s for HD 163296, lower than 
but close to the threshold value noted by lAresu et ali (1201 ll) for 
T Tauri discs. HD 163296 represents an interesting test case for 
the influence of X-rays in Herbig discs. 

We have computed a disc model with parameters identical 
to our preferred model, but with an additional X-ray component 
in the input spectrum. The X-ray luminosity was_ set equal to 
the va lue of Lx — 10 296 erg s 1 observed by iGtinther & Schmittl 
(l2009h for HD 163296. The effect of this on the gas temper- 
ature is illustrated by Fig. [12] (c.f. Fig. [7]). The X-ray heating 
processes produce an extended hot surface gas layer, as noted 
by lAresu et ak - (1201 lh . with gas temperatures >5000K extend- 
ing out to the outer disc. 

There is little discernible effect on the disc chemistry, with 
none of the chemical species masses or predicted line fluxes 
changing by more than a few percent. In all cases the effect of 
the additional X-rays is less than produced by switching to the 
“high UV” input spectrum . We would therefore expect the UV 
to dominate the gas chemistry in a Herbig disc of the sort con- 
sidered here. 


6. Summary and Conclusions 

This paper presents new observations of the far-IR lines of the 
Herbig Ae star HD 163296, obtained using Herschel/PACS as 
part of the GASPS open time key program. These consist of a 
detection of the [Oi] 63. 18 /mi line, as well as upper limits for 
eight additional lines. We have supplemented these observations 
with additional line and continuum data, as well as high resolu- 
tion observations of the variable optical spectra obtained as part 
of the EXPORT collaboration. 

We have computed radiation thermo-chemical disc models 
using the disc code ProDiMo, and employed an evolutionary 
/^-minimisation strategy to find the best simultaneous model fit 
to the continuum and line data. The stellar parameters and UV 
input spectrum for the modelling were determined through de- 
tailed analysis in the UV, optical and near-infrared. 
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Fig. 11. Effect of dust settling on the disc chemical abundances. Left column indicates abundances in a fully mixed model, and 
right column represents a strongly-settled model with otherwise identical parameters. Top row: CO abundance with dust (blue 
dashed lines) and gas (white dashed line) temperature contours. Middle row: LL abundance with 300K gas temperature contour 
(black dashed line) and visual extinction contour (white dashed line). Bottom row: FLO abundance. White dashed lines indicate 
Tgas contours enclosing_^‘hot water” region, red dashed lines indicate contours of UV field strength per hydrogen nucleus, log(y//; H j 
as defined in lDraine & Bertoldii d 1 9961) . enclosing the cool water belt. 
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Fig. 12. Gas temperature structure for the preferred model (col- 
umn 3 in Table [6), irradiated by an additional X-ray spec- 
tral component w i th Lx = 1 0 296 erg s 1 , as observed by 

iGunther & Schmitt! d2009h in HD 163296. This is in compari- 
son to the right hand panel of Fig.|7J in which no X-ray spectral 
component is present. 


We obtain reasonable fits to the observed photometry, ISO- 
SWS spectrum, far-IR line fluxes and millimetre CO profiles for 
a variety of discs, and note that parameter degeneracies preclude 
the precise derivation of the disc properties. In particular, the 
effects of dust settling on the vertical thermal structure of the 
gas in discs strongly influence the line emission, placing further 
limits on the derivation of a disc gas mass from the line data. We 
are able to fit the photometry and line data with a well-mixed 
disc with gas/dust ~ 100, but an equally good fit is possible with 
a disc with gas/dust ~ 22, in which dust settling is present. The 
main advantage of the settled model is that the fit to the observed 
ISO-SWS spectrum is greatly improved, and it seems probable 
that some degree of dust settling and gas-depletion is present 
in the disc. These discs have a radial surface density profile as 
derived by iHughes et al. ( 2008 s). with an exponentially-tapered 
outer edge. 

A model fit to the same data with power-law radial density 
profile gives a gas/dust ratio of ~ 9. This power-law model gives 
the best fit to the observed SED data, although it does very badly 
when it comes to fitting the radial intensity distribution of the 
resolved millimetre continuum. This is due to the flat density 
power index required to fit the SED and the observed CO 3-2 
emission, and presents a challenge for future modelling efforts. 

The emitted line fluxes are in general sensitive to the degree 
of dust settling in the disc. This is a firm result in our models, 
and has serious implications for attempts to derive the disc gas 
mass and other properties from line observations. This settling 
flux enhancement arises from changes to the vertical tempera- 
ture and chemical structure, as settled dust grains allow stellar 
UV to penetrate deeper into the disc. The effect is strongest in 
lines which are formed in the warm gas in the inner disc (e.g. a 
factor ~ 30 increase in the H 2 S(l) line), but the low excitation 
molecular lines are also affected, e.g. a factor ~ 2 increase in the 
13 CO 1-0 line. 


We explore the effects of the observed UV variability in this 
object on the gas chemistry in our models, and conclude that 
this effect is not large enough to affect the observable line fluxes 
beyond the current range of instrumental uncertainty. 

We also examine the effect of X-rays on the gas chemistry 
of our models, and find that while X-rays present a significant 
source of gas-heating in the disc surface layers, the observed X- 
ray luminosity of Lx = 10 29 6 ergs _1 does not significantly alter 
the gas chemistry or line emission. Any effects are smaller than 
those expected as a result of the observed UV variability in this 
object. 

It is difficult to reach any firm conclusions regarding the evo- 
lutionary state of the disc of HD 163296. There is some evidence 
to suggest that the disc is gas-depleted. This is in contrast to the 
result found for the Herbig Ae star HD 169142 by iMeeus et al.l 
( 2010 ). where despite indications that the disc is transitional, it 
was found to be gas-rich. We note that there are uncertainties 
associated with our gas/dust values arising from uncertainties in 
the disc dust composition, and the associated opacity law. This 
is reflected by the factor of ~ 3 spread in the range of derived 


dust masses 

found in the literature for this object (||Natt 

a et al. 

2004; Tann 

rkulam et all 2008a; Mannings & Sargent, 

1997 

Isella et al.. 

2007). All of our derived dust masses are 

within 


the range (5 - 17) x 1 0 4 M 0 from the literature. We also note 
the restrictions placed on our conclusions by the assumption of 
constant dust grain properties throughout the disc, and it has 
been suggested that the d ust properties in discs should in general 
be va riable with radius (Bir nstiel et all 120101 : iGuilloteau et all 

|20TD) . 


It would appear that the modelling of HD 163296 is a far 
from straightforward task, and it is difficult to fit it into the stan- 
dard evolutionary picture developed over the past decade. We 
are unable to find a single disc model which fits perfectly the 
entire wealth of observational data for this object. As well as 
the gas/dust ratio, there is uncertainty regarding the disc flar- 
ing, which itself is strongly tied in to the disc gas heating and 
line emission, as well as being indicative of the disc’s evolution- 
ary state. Evidence of in-falling material and the presence of a 
bipolar outflow in HD 163296 seem to indicate that the star is 
actively accreting material, typical of a young object with a gas- 
rich disc. However, it has been suggested that the observed disc 
geometry and flaring could indicate a star at a later stage of its 
evolution. This would be consistent with the evidence for grain 
growth from our modelling, with all the models requiring large 
grains to fit the observations. There is also possible further evi- 
dence for particle growth from high resolution optical spectra of 
this object (see Appendix). This is also consistent with possible 
evidence for dust settling and gas-depletion from this study, but 
is hard to reconcile with evidence for substantial ongoing activ- 
ity in the inner disc. It is clear that this object is at a fascinating 
stage in its evolution. 


Appendix: echelle spectra analysis 

Sect. 0 provides possible evidence of a low gas-to-dust ratio in 
the HD 163296 circumstellar disc as compared to interstellar 
abundances. This result is also supported in this Appendix from 
the analysis of visible high resolution echelle spectra. 

The EXPORT collaboration (lEiroa et all 1200 ll) obtained 
high resolution echelle spectra (R = A/AA = 48000) in the range 
3800-5900 A during five nights in May and July 1998 using the 
Utrecht Echelle Spectrograph (UES) on the William Herschel 
Telescope, La Palma. 
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May 16 1998 



Radial velocity (km/s) 


Fig. 13. Transient absorption components in 16 May 1998. The 
high resolution observed spectrum of several lines (solid lines) is 
plotted along with the synthetic photospheric spectrum (dashed 
lines; iKuruczl dl993tf ). The short vertical dashes mark the ap- 
proximate location of the metallic-only transient absorptions at 
+50 km/s. A transient red-shifted component of velocity v ~ 
50 km/s is clearly seen in the metallic but not the hydrogen lines. 
In addition a blue-shifted absorption is identified in all lines for 
v 100 km/s. 


Significant spectral variability was found in the optical spec- 
trum. The EXPORT data m onitored those changes on time scales 
of hours, days and months (iMendigutla et alll201 ll) . Almost any 
strong photospheric line is variable to some extent. In particular. 
Transient Absorption Components (TACs) are observed in many 
hydrogen (Balmer series) and metallic lines, such as Ca n K & 
H, Nai D2 & D1 and the Fen 4924, 5018 and 5169 A triplet. 
The Appendix summarizes the main results on the dynamics of 
the circumstellar gas obtained from the TACS. 

Many TACs obse rved are similar to tho se studied by 
iNatta et afl (l2000h and iMora et all d2002l 120041) in similar age 
Herbig Ae stars. These authors propose magnetospheric accre- 
tion of primordial gas as the most likely origin. That is, the TACs 
in Herbig Ae stars typically trace the kinematics of material with 
gas to dust ratios similar to that of the interstellar medium. 

However, a distinctive feature of HD 163296 is the presence 
of low velocity (v ~ 50 km/s) red-shifted TACs only observed 
in refractory elements: Fei, Feii, Tin and Sen. These “metal- 
lic” TACs have no clear counterpart neither in H nor Ca n and 
Nai. They were observed during all nights except 28 July 1998. 


July 16 1998 



Radial velocity (km/s) 


Fig. 14. Transient absorption components in 16 July 1998. The 
high resolution observed spectrum of several lines (solid lines) 
is plotted along with the photospheric spectrum (dashed lines). 
No clear transient absorption is detected in the metallic lines. On 
the other hand a significant blue-shifted absorption is detected in 
all Balmer lines at v ~ -100 km/s. 


These transient absorption s could be better expl a ined i n terms 
of the model developed by lLagrange-Henri et al.l ( 1 988 ) for the 
evaporation of solid b odies in the H Pic d e bris d isc and subse- 
quently postulated by iGrinin et al.l (1199 ll 1 1994b for accretion 
discs around Herbig Ae/Be stars. 

Fig. fl3l shows selected line profiles observed in 16 May 
1998. Wavelengths have been converted to radial velocities in 
the stellar reference system. The approximate location of the 
metallic-only transient absorptions at +50 km/s is displayed 
with a short vertical dash. The underlying photosphere is rep- 
resented by a synthetic Kurucz model computed with T e g = 
9250 K, log g* = 4.07 and [Fe/H]=+0.20, using the ATLAS9 
and SYNTHE codes by IKuruczl (Il993h . In this work we have 
used the GNU Linux version of the codes available onlintfl 
(ISbordone et all l2004t) . It is clear that no Balmer line absorp- 
tion counterpart can be identified, even for the higher Hd and 
HC lines in the series for which no significant core line emis- 
sion is apparent. On the other hand, blue-shifted absorption with 
v ~ -100 km/s is clearly seen both in metallic and hydrogen 
lines. The metallic only transient events are not an artifact gener- 
ated by poor photosphere characterisation, as revealed by Fig. M 


http://wwwuser.oat.ts.astro.it/atmos/ 
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(electronic edition only), which shows the spectra obtained in 28 
July 1998, the only night where no metallic only events where 
identified. 

Similar results were found for WW Vul by iMora et aD 
d2004j> . but not for the other four stars studied. In addition, both 
hydrogen-rich and hydrogen-less TACs are simultaneously de- 
tected in HD 163296 and WW Vul. If particle growth indeed 
explains the metallic-only events, that would suggest the circum- 
stellar disc around HD 163296 (and incidentally WW Vul) has 
an evolutionary status between a pure primordial gas rich accre- 
tion disc and a gas-depleted second generation debris disc. That 
finding is fully consistent with the results obtained in Sec.Ofrom 
the analysis of Herschel/PACS far-IR spectroscopy and ground- 
based mm CO observations. 

Kinematics of the circumstellar gas 

A detailed study of the dynamics of the circumstellar gas around 
HD 163296 using the method presented in iMora et al.l d2004l) 
is out the scope of this paper. The main results are, however, 
summarized here. 

The following lines were studied: the Balmer series (H/J, Hy, 
H<5, He and He), The Nai D and Can HK doublets, the Fen 
4924, 5018 and 5169 A triplet, Fe i 4046 A, Ti n 4444 and 4572 
A and Sc ii 4247 A. Three types of variability were always si- 
multaneously observed: core line emission (most prominent in 
the Balmer, Na i and Ca n lines), transient red-shifted absorptions 
(associated to in-falling material, observed in all lines) and tran- 
sient blue-shifted absorptions (associated to outflows, observed 
in all lines). HD 163296 wa s more active than the stars studied 
by IMora et all (I2002L |2004|) . In particular, the larger line core 
emission made the analysis more difficult. 

The in-falling material always span a large range of veloci- 
ties, from near zero to +300 km/s for the red end of the higher 
velocity components. The stellar winds maximum velocities can 
also be high. A maximum of -400 km/s was observed in 17 May 
1998. 

Fig. IT51 shows a summary of all the transient absorptions 
identified. Features with similar radial velocity observed in dif- 
ferent lines are assumed to share a common origin. The average 
radial velocity and the standard deviation are displayed as points 
and error bars in the figure. The number of lines used to com- 
pute the average is given next to each point. The fainter lines are 
given half weight in the average (see lMora et alJ ( 2004 1 for addi- 
tional details). Points with similar velocity on nearby nights are 
also assumed to be produced by the same clump of circumstel- 
lar material and are grouped in ’’events”. Events are displayed 
in the figure with boxes enclosing a certain number of average 
points. Red boxes represent in-falling material, blue boxes out- 
flows and green boxes in-falling material only detected in metal- 
lic lines. Metallic-only events have smaller radial velocities than 
hydrogen-rich in-falling material. 

One of the most interesting events observed is #4, a strong 
wind observed in plenty of Balmer and metallic lines during July 
1998. The most relevant feature is an almost constant accelera- 
tion rate of ~ -0.4 m/s 2 . The average radial velocity starts at -120 
km/s on 28 July 1998 and grows up to -235 km/s in 31 July 
1998. The width of the absorption features has an almost con- 
stant value of 120 km/s, supporting the hypothesis of a common 
origin of all the transient absorptions in the event. 
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Fig. 15. Kinematics of the circumstellar gas. Features with similar radial velocity observed in different lines and different dates are 
assumed to share a common origin and grouped into events. Each event is represented by a box enclosing a number of points. Each 
point represents the weighted average of the radial velocity simultaneously measured for several lines. The error bar is the standard 
deviation. The weighted number of lines is given next to each point. Red boxes represent in-falling material, blue boxes out-flowing 
material and green boxes in-falling material only detected in metallic lines. 
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